2 Read Bismark files

2.1 Negative control

## Reading ./input/NC_rep1.cov.gz 
## Reading ./input/NC_rep2.cov.gz 
## Reading ./input/NC_rep3.cov.gz 
## Reading ./input/NC_rep4.cov.gz 
## Reading ./input/NC_rep5.cov.gz 
## Reading ./input/NC_rep6.cov.gz 
## Hashing ...
## Collating counts ...

2.2 Simulated data

## Reading ./input/sim_rep1.cov.gz 
## Reading ./input/sim_rep2.cov.gz 
## Reading ./input/sim_rep3.cov.gz 
## Reading ./input/sim_rep4.cov.gz 
## Reading ./input/sim_rep5.cov.gz 
## Reading ./input/sim_rep6.cov.gz 
## Hashing ...
## Collating counts ...

3 Filtering and plotting

3.1 Function

beta_calc <- function(objectDGE, coverageThreshold, anno = anno) {
  ## Function to calculate Beta values from DGE object obtained from readBismark2DGE function.
  ## objectDGE: DGEList
  ## coverageThreshold: Coverage threshold to be used for all samples
  ## Return: Coverage, Coverage melted table, distribution, Beta value table, Beta value melted table, distribution

  # Counts (raw coverage)
  df <- data.frame(objectDGE$counts, stringsAsFactors = F, check.names = F)

  me.raw <- df[, grep("Me", colnames(df))]
  un.raw <- df[, grep("Un", colnames(df))]
  cov.raw <- me.raw + un.raw


  # Filter based on coverage
  data <- NULL

  if (is.null(coverageThreshold)) {
    data <- df
    coverageThreshold <- "NO filtering"
  } else {
    data <- df[rowSums(cov.raw > coverageThreshold) == ncol(cov.raw), ]
  }

  # Beta valus after filtering
  me <- data[, grep("Me", colnames(data))]
  un <- data[, grep("Un", colnames(data))]
  beta <- me / (me + un)
  colnames(beta) <- as.character(sapply(colnames(beta), function(x) strsplit(x, "-")[[1]][1]))

  # Beta values for plotting
  reshape.beta <- data.frame(t(beta))
  reshape.beta$Samples <- rownames(reshape.beta)
  reshape.beta.melt <- melt(reshape.beta)

  # Plotting Beta value
  p1 <- ggplotly(ggplot(reshape.beta.melt, aes(x = value)) +
    geom_histogram(aes(y = ..density..), alpha = 0.7, fill = "#333333") +
    geom_density(fill = "#CC6666", alpha = 0.3) +
    theme(panel.background = element_rect(fill = "#ffffff")) +
    facet_wrap(~Samples, nrow = 3, shrink = T) +
    ggtitle(paste0("Beta distribution (coverage > ", coverageThreshold, ")")) +
    labs(x = "Beta value", y = ""))

  # Coverage plot after filtering
  cov <- me + un
  colnames(cov) <- as.character(sapply(colnames(cov), function(x) strsplit(x, "-")[[1]][1]))
  reshape.cov <- data.frame(t(cov))
  reshape.cov$Samples <- rownames(reshape.cov)
  reshape.cov.melt <- melt(reshape.cov)
  reshape.cov.melt.log2 <- reshape.cov.melt
  reshape.cov.melt.log2$value <- log2(reshape.cov.melt.log2$value)

  p2 <- ggplotly(ggplot(reshape.cov.melt.log2, aes(x = value)) +
    geom_histogram(aes(y = ..density..), alpha = 0.7, fill = "#333333") +
    geom_density(fill = "#20AAEA", alpha = 0.3) +
    theme(panel.background = element_rect(fill = "#ffffff")) +
    facet_wrap(~Samples, nrow = 3, shrink = T) +
    ggtitle(paste0("Coverage distribution (coverage > ", coverageThreshold, ")")) +
    labs(x = "Coverage (log2 scaled)", y = ""))

  # PCA plot

  ## Coverage
  pca.cov <- prcomp(t(cov))

  p3.cov <- plot_ly(
    x = pca.cov$x[, 1], y = pca.cov$x[, 2],
    data = anno, color = ~Group, mode = "markers",
    symbol = ~ as.numeric(as.factor(Group)),
    marker = list(size = 10), hoverinfo = "text",
    text = ~ paste("<b>Group:</b> ", Group, "<br><b>SampleID:</b> ", names),
    type = "scatter"
  ) %>%
    layout(
      title = paste0(
        "<b>Coverage PCA plot (<span style='color: blue;'>coverage > ",
        coverageThreshold, "</span>)<b>"
      ),
      xaxis = list(title = "<b><i>PC1<i><b>"),
      yaxis = list(title = "<b><i>PC2<i><b>")
    )


  ## Beta
  pca.beta <- prcomp(t(beta))

  p3.beta <- plot_ly(
    x = pca.beta$x[, 1], y = pca.beta$x[, 2],
    data = anno, color = ~Group, mode = "markers",
    symbol = ~ as.numeric(as.factor(Group)),
    marker = list(size = 10), hoverinfo = "text",
    text = ~ paste("<b>Group:</b> ", Group, "<br><b>SampleID:</b> ", names),
    type = "scatter"
  ) %>%
    layout(
      title = paste0(
        "<b>Beta PCA plot (<span style='color: blue;'>coverage > ",
        coverageThreshold, "</span>)<b>"
      ),
      xaxis = list(title = "<b><i>PC1<i><b>"),
      yaxis = list(title = "<b><i>PC2<i><b>")
    )


  # MDS plot

  ## Coverage
  mds.cov <- cmdscale(dist(t(cov)))

  p4.cov <- plot_ly(
    x = mds.cov[, 1], y = mds.cov[, 2],
    data = anno, color = ~Group, mode = "markers",
    symbol = ~ as.numeric(as.factor(Group)),
    marker = list(size = 10), hoverinfo = "text",
    text = ~ paste("<b>Group:</b> ", Group, "<br><b>SampleID:</b> ", names),
    type = "scatter"
  ) %>%
    layout(
      title = paste0(
        "<b>Coverage MDS plot (<span style='color: blue;'>coverage > ",
        coverageThreshold, "</span>)<b>"
      ),
      xaxis = list(title = "<b><i>MDS1<i><b>"),
      yaxis = list(title = "<b><i>MDS2<i><b>")
    )


  ## Beta
  mds.beta <- cmdscale(dist(t(beta)))

  p4.beta <- plot_ly(
    x = mds.beta[, 1], y = mds.beta[, 2],
    data = anno, color = ~Group, mode = "markers",
    symbol = ~ as.numeric(as.factor(Group)),
    marker = list(size = 10), hoverinfo = "text",
    text = ~ paste("<b>Group:</b> ", Group, "<br><b>SampleID:</b> ", names),
    type = "scatter"
  ) %>%
    layout(
      title = paste0(
        "<b>Beta MDS plot (<span style='color: blue;'>coverage > ",
        coverageThreshold, "</span>)<b>"
      ),
      xaxis = list(title = "<b><i>MDS1<i><b>"),
      yaxis = list(title = "<b><i>MDS2<i><b>")
    )

  # SVD

  ## Coverage

  svd.cov <- svd(x = t(cov))

  p5.cov <- plot_ly(
    x = svd.cov$u[, 1], y = svd.cov$u[, 2],
    data = anno, color = ~Group, mode = "markers",
    symbol = ~ as.numeric(as.factor(Group)),
    marker = list(size = 10), hoverinfo = "text",
    text = ~ paste("<b>Group:</b> ", Group, "<br><b>SampleID:</b> ", names),
    type = "scatter"
  ) %>%
    layout(
      title = paste0(
        "<b>Coverage SVD plot (<span style='color: blue;'>coverage > ",
        coverageThreshold, "</span>)<b>"
      ),
      xaxis = list(title = "<b><i>SVD1<i><b>"),
      yaxis = list(title = "<b><i>SVD2<i><b>")
    )

  ## Beta
  svd.beta <- svd(x = t(beta))


  p5.beta <- plot_ly(
    x = svd.beta$u[, 1], y = svd.beta$u[, 2],
    data = anno, color = ~Group, mode = "markers",
    symbol = ~ as.numeric(as.factor(Group)),
    marker = list(size = 10), hoverinfo = "text",
    text = ~ paste("<b>Group:</b> ", Group, "<br><b>SampleID:</b> ", names),
    type = "scatter"
  ) %>%
    layout(
      title = paste0(
        "<b>Beta SVD plot (<span style='color: blue;'>coverage > ",
        coverageThreshold, "</span>)<b>"
      ),
      xaxis = list(title = "<b><i>SVD1<i><b>"),
      yaxis = list(title = "<b><i>SVD2<i><b>")
    )


  # Clustering

  ## Coverage
  p6.cov <- ggdendrogram(hclust(dist(t(cov)))) +
    ggtitle(paste0("Coverage dendrogram (coverage > ", coverageThreshold, ")"))

  ## Beta
  p6.beta <- ggdendrogram(hclust(dist(t(beta)))) +
    ggtitle(paste0("Beta dendrogram (coverage > ", coverageThreshold, ")"))


  list <- list(
    cov, reshape.cov.melt, p2, beta, reshape.beta.melt, p1,
    p3.cov, p4.cov, p5.cov, p6.cov, p3.beta, p4.beta, p5.beta, p6.beta
  )
  names(list) <- c(
    "Coverage", "reshapedCoverage", "DHPlotCov", "BetaValues", "reshapedBetaValues", "DHplotBeta",
    "PCAcov", "MDScov", "SVDcov", "HCcov", "PCAbeta", "MDSbeta", "SVDbeta", "HCbeta"
  )

  return(list)
}

4 Coverage Plots

4.1 Negative control

4.2 Simulated data

5 Beta Plots

5.1 Negative control

5.2 Simulated data

6 Clustering of samples

6.1 Negative control

6.2 Simulated data

7 PCA Plots

7.1 Negative control

7.2 Simulated data

8 MDS Plots

8.1 Negative control

8.2 Simulated data

9 SVD Plots

9.1 Negative control

9.2 Simulated data

10 SessionInfo

## ─ Session info ──────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.5.1 (2018-07-02)
##  os       Ubuntu 16.04.5 LTS          
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       Europe/Zurich               
##  date     2019-01-08                  
## 
## ─ Packages ──────────────────────────────────────────────────────────────
##  package      * version  date       lib source        
##  abind          1.4-5    2016-07-21 [1] CRAN (R 3.5.1)
##  ade4           1.7-13   2018-08-31 [1] CRAN (R 3.5.1)
##  ape            5.2      2018-09-24 [1] CRAN (R 3.5.1)
##  assertthat     0.2.0    2017-04-11 [1] CRAN (R 3.5.1)
##  autoplotly   * 0.1.2    2018-04-21 [1] CRAN (R 3.5.1)
##  backports      1.1.3    2018-12-14 [1] CRAN (R 3.5.1)
##  bindr          0.1.1    2018-03-13 [1] CRAN (R 3.5.1)
##  bindrcpp     * 0.2.2    2018-03-29 [1] CRAN (R 3.5.1)
##  Cairo          1.5-9    2015-09-26 [1] CRAN (R 3.5.1)
##  callr          3.1.1    2018-12-21 [1] CRAN (R 3.5.1)
##  cli            1.0.1    2018-09-25 [1] CRAN (R 3.5.1)
##  cluster        2.0.7-1  2018-04-09 [1] CRAN (R 3.5.1)
##  ClusterR     * 1.1.7    2018-12-10 [1] CRAN (R 3.5.1)
##  colorspace     1.3-2    2016-12-14 [1] CRAN (R 3.5.1)
##  crayon         1.3.4    2017-09-16 [1] CRAN (R 3.5.1)
##  crosstalk      1.0.0    2016-12-21 [1] CRAN (R 3.5.1)
##  data.table   * 1.11.8   2018-09-30 [1] CRAN (R 3.5.1)
##  desc           1.2.0    2018-05-01 [1] CRAN (R 3.5.1)
##  devtools       2.0.1    2018-10-26 [1] CRAN (R 3.5.1)
##  digest         0.6.18   2018-10-10 [1] CRAN (R 3.5.1)
##  dplyr          0.7.8    2018-11-10 [1] CRAN (R 3.5.1)
##  edgeR        * 3.22.5   2018-10-02 [1] Bioconductor  
##  evaluate       0.12     2018-10-09 [1] CRAN (R 3.5.1)
##  FD             1.0-12   2014-08-19 [1] CRAN (R 3.5.1)
##  fs             1.2.6    2018-08-23 [1] CRAN (R 3.5.1)
##  geometry       0.3-6    2015-09-09 [1] CRAN (R 3.5.1)
##  ggdendro     * 0.1-20   2016-04-27 [1] CRAN (R 3.5.1)
##  ggfortify    * 0.4.5    2018-05-26 [1] CRAN (R 3.5.1)
##  ggplot2      * 3.1.0    2018-10-25 [1] CRAN (R 3.5.1)
##  glue           1.3.0    2018-07-17 [1] CRAN (R 3.5.1)
##  gmp            0.5-13.2 2018-07-14 [1] CRAN (R 3.5.1)
##  gridExtra      2.3      2017-09-09 [1] CRAN (R 3.5.1)
##  gtable         0.2.0    2016-02-26 [1] CRAN (R 3.5.1)
##  gtools       * 3.8.1    2018-06-26 [1] CRAN (R 3.5.1)
##  hms            0.4.2    2018-03-10 [1] CRAN (R 3.5.1)
##  htmltools      0.3.6    2017-04-28 [1] CRAN (R 3.5.1)
##  htmlwidgets    1.3      2018-09-30 [1] CRAN (R 3.5.1)
##  httpuv         1.4.5.1  2018-12-18 [1] CRAN (R 3.5.1)
##  httr           1.4.0    2018-12-11 [1] CRAN (R 3.5.1)
##  jpeg           0.1-8    2014-01-23 [1] CRAN (R 3.5.1)
##  jsonlite       1.6      2018-12-07 [1] CRAN (R 3.5.1)
##  knitr          1.21     2018-12-10 [1] CRAN (R 3.5.1)
##  labeling       0.3      2014-08-23 [1] CRAN (R 3.5.1)
##  later          0.7.5    2018-09-18 [1] CRAN (R 3.5.1)
##  lattice        0.20-38  2018-11-04 [1] CRAN (R 3.5.1)
##  lazyeval       0.2.1    2017-10-29 [1] CRAN (R 3.5.1)
##  limma        * 3.36.5   2018-09-20 [1] Bioconductor  
##  locfit         1.5-9.1  2013-04-20 [1] CRAN (R 3.5.1)
##  magic          1.5-9    2018-09-17 [1] CRAN (R 3.5.1)
##  magrittr       1.5      2014-11-22 [1] CRAN (R 3.5.1)
##  MASS           7.3-51.1 2018-11-01 [1] CRAN (R 3.5.1)
##  Matrix         1.2-15   2018-11-01 [1] CRAN (R 3.5.1)
##  memoise        1.1.0    2017-04-21 [1] CRAN (R 3.5.1)
##  mgcv           1.8-26   2018-11-21 [1] CRAN (R 3.5.1)
##  mime           0.6      2018-10-05 [1] CRAN (R 3.5.1)
##  munsell        0.5.0    2018-06-12 [1] CRAN (R 3.5.1)
##  nlme           3.1-137  2018-04-07 [1] CRAN (R 3.5.1)
##  OpenImageR     1.1.3    2018-12-09 [1] CRAN (R 3.5.1)
##  permute        0.9-4    2016-09-09 [1] CRAN (R 3.5.1)
##  pillar         1.3.1    2018-12-15 [1] CRAN (R 3.5.1)
##  pkgbuild       1.0.2    2018-10-16 [1] CRAN (R 3.5.1)
##  pkgconfig      2.0.2    2018-08-16 [1] CRAN (R 3.5.1)
##  pkgload        1.0.2    2018-10-29 [1] CRAN (R 3.5.1)
##  plotly       * 4.8.0    2018-07-20 [1] CRAN (R 3.5.1)
##  plyr           1.8.4    2016-06-08 [1] CRAN (R 3.5.1)
##  png            0.1-7    2013-12-03 [1] CRAN (R 3.5.1)
##  prettyunits    1.0.2    2015-07-13 [1] CRAN (R 3.5.1)
##  processx       3.2.1    2018-12-05 [1] CRAN (R 3.5.1)
##  promises       1.0.1    2018-04-13 [1] CRAN (R 3.5.1)
##  ps             1.3.0    2018-12-21 [1] CRAN (R 3.5.1)
##  purrr          0.2.5    2018-05-29 [1] CRAN (R 3.5.1)
##  R6             2.3.0    2018-10-04 [1] CRAN (R 3.5.1)
##  RColorBrewer   1.1-2    2014-12-07 [1] CRAN (R 3.5.1)
##  Rcpp           1.0.0    2018-11-07 [1] CRAN (R 3.5.1)
##  readr          1.3.1    2018-12-21 [1] CRAN (R 3.5.1)
##  remotes        2.0.2    2018-10-30 [1] CRAN (R 3.5.1)
##  reshape2     * 1.4.3    2017-12-11 [1] CRAN (R 3.5.1)
##  rlang          0.3.0.1  2018-10-25 [1] CRAN (R 3.5.1)
##  rmarkdown      1.11     2018-12-08 [1] CRAN (R 3.5.1)
##  rprojroot      1.3-2    2018-01-03 [1] CRAN (R 3.5.1)
##  scales         1.0.0    2018-08-09 [1] CRAN (R 3.5.1)
##  sessioninfo    1.1.1    2018-11-05 [1] CRAN (R 3.5.1)
##  shiny          1.2.0    2018-11-02 [1] CRAN (R 3.5.1)
##  stringi        1.2.4    2018-07-20 [1] CRAN (R 3.5.1)
##  stringr        1.3.1    2018-05-10 [1] CRAN (R 3.5.1)
##  testthat       2.0.1    2018-10-13 [1] CRAN (R 3.5.1)
##  tibble         1.4.2    2018-01-22 [1] CRAN (R 3.5.1)
##  tidyr          0.8.2    2018-10-28 [1] CRAN (R 3.5.1)
##  tidyselect     0.2.5    2018-10-11 [1] CRAN (R 3.5.1)
##  tiff           0.1-5    2013-09-04 [1] CRAN (R 3.5.1)
##  usethis        1.4.0    2018-08-14 [1] CRAN (R 3.5.1)
##  vegan          2.5-3    2018-10-25 [1] CRAN (R 3.5.1)
##  viridis      * 0.5.1    2018-03-29 [1] CRAN (R 3.5.1)
##  viridisLite  * 0.3.0    2018-02-01 [1] CRAN (R 3.5.1)
##  withr          2.1.2    2018-03-15 [1] CRAN (R 3.5.1)
##  xfun           0.4      2018-10-23 [1] CRAN (R 3.5.1)
##  xtable         1.8-3    2018-08-29 [1] CRAN (R 3.5.1)
##  yaml           2.2.0    2018-07-25 [1] CRAN (R 3.5.1)
## 
## [1] /home/ubuntu/R/x86_64-pc-linux-gnu-library/3.5
## [2] /usr/local/lib/R/site-library
## [3] /usr/lib/R/site-library
## [4] /usr/lib/R/library